Last year, the number of software developers was estimated to be 26.4 million. The field thrives on open source tools & technologies, with new tools created every day. Given the sheer number of tools available, it's very possible that millions of unique combinations of these tools and technologies exist amongst software developers.¶
This project applies Data Cleaning, Exploratory Data Analysis & Machine Learning techniques to better understand the connection between a developer's tools and their profile type.¶
The Presentation¶
Overview of the Data
- Top Skills
- Top Profiles
Defining our Dataset
- Database
- Machine Learning (ML)
- Quality Assurance (QA)
- Frontend
- Fullstack
Transforming the Data
Reducing Noise with Feature Selection
The Machine Learning Models¶
Logistic Regression¶
Class Balance
Benchmark Model
Scoring with the ROC Curve
Class Imbalance with PySpark
Model Performance Metrics
Binary Decision Tree¶
Benchmark Model
Model Performance with the Classification Report
Scoring with ROC_AUC
Visualizing the Decision Tree
Model Performance Metrics
Multinomial (Multi-class) Random Forest¶
Benchmark Model
Model Performance with the Classification ReportThe Data¶
data_overview = pd.read_csv("../Celerative_Working/Data_Cleaning/data_overview.csv")
data_overview.drop(["Unnamed: 0"], axis=1, inplace=True)
data_overview.head(5)
How An Algorithm sees the Data¶
| Observations or Samples | Classes, Predicted vs. Actual (y) | Features (X) |
|---|---|---|
| observation #1 | profile_A | skill_1, skill_5, skill_3... |
| observation #2 | profile_B | skill_5, skill_2, skill_4... |
| observation #3 | profile_C | skill_4, skill1_2, skill_1... |
| observation #4 | profile_B | skill_4, skill1_1, skill_5... |
| observation #5 | profile_A | skill_3, skill1_4, skill_5... |
Data Cleaning¶
Noise is data which does not add meaningful information for the algorithm.¶
Controlling for noise yields more reliable predictive analysis when applying Machine Learning algorithms.¶
This is an example of noisy data. Left unchecked, this could cause the algorithm to return results which don't acurately reflect the true data:¶
| Similar | Strings | |
|---|---|---|
| 'apache' | 'apache spark' | |
| 'apache ant' | 'apache tomcat' | |
| 'rest' | 'rest apis' | |
| 'android studio' | 'android' | 'android sdk' |
| 'asp.net mvc' | 'asp.net' | 'asp' |
| 'visual studio code' | 'microsoft visual studio' | |
| 'angular js' | 'angular' | |
| 'c#' | 'c' | 'c++' |
| 'asp.net mvc' | 'asp.net' | |
| 'google app engine' | 'google cloud' | 'google compute engine' |
| 'google analytics' | 'google adsense' | 'google adwords', |
| 'phpstorm' | 'php' | |
| 'net' | 'net framework' | |
| 'aws ec2' | 'aws s3 |
Before¶
print("Before the data was cleaned and processed, the dataset consisted of:\n")
print(
"\t - ",
data_overview.profile.nunique(),
'unique profile types or "classes" \n',
)
print("\t & ")
print(
"\n\t - ",
data_overview.skills.explode().nunique(),
'skills or "features"',
)
After¶
print(
"Using Natural Language Processing and Text Pre-Processing techniques provided via the Scikit-Learn library, skills & profile strings were standardized to:"
)
print("\t - ", clean_talentpool.profile.nunique(), "profiles (classes)")
print("\t & ")
print("\t - ", clean_talentpool.skills.explode().nunique(), "skills (features)")
Top Skills¶
# Generating wordcloud
wordcloud = WordCloud(
width=3000,
height=2000,
random_state=1,
max_words=250,
background_color="black",
colormap="Blues",
collocations=False,
stopwords=STOPWORDS,
).generate(
series_to_str(clean_talentpool["skills"])
) # Plot
plot_cloud(wordcloud)
Top Profiles¶
pd.DataFrame(clean_talentpool["profile"].value_counts())[1:19]
It's clear that a few of the top profile types are described by their key skills only, such as "frontend php javascript" , "javascript frontend design", etc. These top "keywords" may reveal more than about the data's profiles when isolated from non-descritptive words like "developer" and "engineer".
# Drop empty row & subsetting for top 25 keywords
dotplot.drop(14, inplace=True)
dotplot.reset_index(drop=True)[:25]
# Draw plot
y = dotplot.profile_keywords[:45]
fig, ax = plt.subplots(figsize=(10, 12), dpi=80)
ax.hlines(
y=y,
xmin=10,
xmax=205,
color="gray",
alpha=0.7,
linewidth=1,
linestyles="dashdot",
)
ax.scatter(
y=y,
x=dotplot.loc[:45, ["count"]],
s=75,
color="firebrick",
alpha=0.7,
)
# Title, Label, Ticks and Ylim
ax.set_title("Most Common Profile Keywords", fontdict={"size": "22"})
ax.set_xlabel("Count of Keywords", fontdict={"size": "22"})
ax.set_yticklabels(y, fontdict={"horizontalalignment": "right", "size": "14"})
ax.set_xlim(5, 200)
plt.show()
Defining our Dataset¶
For the purposes of this presentation, we will narrow our focus to look at machine learning model performance using a dataset made up of five (5) profiles:¶
- "Database"
- "ML" combined those profiles containing the phrase "machine learning"
- "QA" combined those profiles containing the phrase "quality assurance"
- "Fullstack"
- "Frontend"
Makeup of the Dataset¶
# Get percentage of developer profiles
(five_profiles.profile.value_counts(normalize=True) * 100).astype(int)
print(
"Although profile types have been cut to 5 broad profile categories, \
\nthere are still",
ex_skills.skills.nunique(),
"skills (vs. the 812 seen with all profiles). \
\nSo although many of the classes have been cut, most skills remain,\
\nqualifying the exploration of this particular subset of the data as valid.",
)
Transforming the Data¶
Our Data's Format¶
| Observations/Samples | Classes (y) | Features (X) |
|---|---|---|
| observation #1 | profile_A | skill_1, skill_5, skill_3... |
| observation #2 | profile_B | skill_5, skill_2, skill_4... |
| observation #3 | profile_C | skill_4, skill1_2, skill_1... |
| observation #4 | profile_B | skill_4, skill1_1, skill_5... |
| observation #5 | profile_A | skill_3, skill1_4, skill_5... |
Machine Learning Data Format¶
| Observations/Samples | Classes (y) | Feature1 (X) | Feature2 (X) | Feature3 (X) | Feature4 (X)... |
|---|---|---|---|---|---|
| observation #1 | profile_A | skill_1 | skill_2 | skill_3 | skill_4 ... |
| observation #2 | profile_B | skill_1 | skill_2 | skill_3 | skill_4 ... |
| observation #3 | profile_C | skill_1 | skill_2 | skill_3 | skill_4 ... |
| observation #4 | profile_B | skill_1 | skill_2 | skill_3 | skill_4 ... |
| observation #5 | profile_A | skill_1 | skill_2 | skill_3 | skill_4 ... |
The Sci-Kit Learn MultiLabelBinarizer() Object¶
sklearn's MultiLabelBinarizer() Object is used to take each skill and expand it, so that each skill (feature) has its own column. Should that sample, or developer, have the specified skill, the corresponding row is marked with a 1. Otherwise, it's marked with a 0.
This effectively converts all skills to binary variables, which allows the features (skills) to be represented on the x-axis, while observations are represented on the y-axis. Machine Learning algorithms need data formatted in this way in order to process it correctly.
The Result¶
df = binarized(ex_skills)
df.iloc[:, 1:8]
X¶
is recognized as the variable which stores the "input variables", or features on which the algorithm is trained.
Here features are assigned to X:
# Assign features to X
X = df.iloc[:, 5:]
X
y¶
is recognized as the variable which stores class labels in classification problems.
Here class labels are assigned to y:
# Assign target class to y
y = df.profile
y
Now that the data has been formatted to be understandable to the algorithm, inconsistencies or noise in the data should be addressed.¶
Feature Engineering or Feature Selection¶
involves choosing features of the data that are more relevant to the model performing well.
This makes data input to the model less noisy, lets the model train faster, simplifies the model for interpretation, and improves generalization to outside data.
- Recursive Feature Elimination (RFE)¶
- Correlation¶
Correlation¶
# Create the correlation matrix
corr_matrix = X.corr()
# Subset the matrix
corr_matrix.iloc[:11, :11]
# Selecting for those columns with correlation relationship less than 50%
columns = np.full((corr_matrix.shape[0],), True, dtype=bool)
for i in range(corr_matrix.shape[0]):
for j in range(i + 1, corr_matrix.shape[0]):
if corr_matrix.iloc[i, j] >= 0.5:
if columns[j]:
columns[j] = False
By selecting for features more than 50% dissimilar, features have been reduced from 526 to 278, effectively reducing the "dimensionality" or column vectors of the dataset. This is one of an array of techniques used to select for the most important features of a dataset.
# newly selected features
X
Recursive Feature Elimination (RFE)¶
Dimensionality Reduction, Matrix style¶
# Performing RFE
rForest = RandomForestClassifier()
rfe = RFE(rForest).fit(X, y)
# Rank of 1 indicates the features selected
rfe.ranking_
from IPython.display import IFrame
IFrame(
"https://giphy.com/embed/AOSwwqVjNZlDO",
width="100%",
height="600",
scrolling="no",
frameborder="0",
)
Optimizing Recursive Feature Elimination (RFE):
The number of features returning the best balanced_accuracy score is 130 features. This is the optimal number of features at which the random forest algorithm performs best! 
RFE chooses features of the data that are more relevant to the model performing well. Doing so makes the data input to the model less noisy, lets the model train faster, simplifies the model for interpretation, and improves generalisation to outside data.
Use of K_Means methods to find patterns in the data¶
Use of the elbow method reveals the optimal number of cluster by which to cluster to data:¶
from yellowbrick.cluster import KElbowVisualizer
# Instantiate the clustering model and visualizer
model = KMeans()
visualizer = KElbowVisualizer(model, k=(1, 130))
visualizer.fit(X_T) # Fit the data to the visualizer
visualizer.show()
DBSCAN (Density-based spatial clustering of applications with noise)¶
is actually an unsupervised machine learning algorithm by itself, but is often used as an Exploratory Data Analysis technique as it reveals natural clusters in the data. The graph above has grouped the remaining 130 features (or skills) into 16 clusters. Those will be used below to find features of importance.
- DBSCAN does not require us that the number of clusters be specified a-priori (beforehand), avoids outliers, and works quite well with arbitrarily shaped and sized clusters.
The silhouette ranges from −1 to +1, with a high value indicating that the object is well matched to its own cluster and poorly matched to neighboring clusters. So a value closer to 0, like ours here, indicates the clusters aren't well separated and share similarity to other clusters, which matches what we see with feature variance in that not many features are causing a high amount of variance in the predictive analysis.
kmeans = KMeans(n_clusters=34).fit(X_130)
normalized_vectors = preprocessing.normalize(X_130)
normalized_kmeans = KMeans(n_clusters=34).fit(normalized_vectors)
min_samples = X_130.shape[1] + 1
dbscan = DBSCAN(eps=3.5, min_samples=min_samples).fit(X_130)
# Print results
print(
"kmeans: {0:10.4f}".format(
silhouette_score(X_130, kmeans.labels_, metric="euclidean")
)
)
print(
"Cosine kmeans:{0:10.4f}".format(
silhouette_score(
normalized_vectors, normalized_kmeans.labels_, metric="cosine"
),
)
)
The DBSCAN algorithm should be used to find associations and structures in data that are hard to find manually but that can be relevant and useful to find patterns and predict trends. Below, we select the top ten (10) variables with the highest variance. A high variance indicates that this variable is likely to play an important role in creating the clusters.
results = pd.DataFrame(columns=["Variable", "Var"])
for column in df_mean.columns[1:]:
results.loc[len(results), :] = [column, np.var(df_mean[column])]
selected_columns = (
list(
results.sort_values(
"Var",
ascending=False,
)
.head(10)
.Variable.values
)
+ ["dbscan"]
)
tidy = df_scaled[selected_columns].melt(id_vars="dbscan")
sns.barplot(x="dbscan", y="value", hue="variable", data=tidy)
print(X_130.columns.to_numpy()[66])
print(X_130.columns.to_numpy()[89])
print(X_130.columns.to_numpy()[90])
print(X_130.columns.to_numpy()[91])
print(X_130.columns.to_numpy()[92])
print(X_130.columns.to_numpy()[93])
print(X_130.columns.to_numpy()[94])
print(X_130.columns.to_numpy()[95])
print(X_130.columns.to_numpy()[96])
Combining Clustering & Correlation techniques to visualize patterns amongst 75 features¶
# Set y_ticks
y_axis_labels = list(df.index)
# Normalize the data
df_norm_col = (df - df.mean()) / df.std()
# Plot the data
ax = sns.clustermap(df_norm_col, yticklabels=y_axis_labels)
Factor Analysis¶
X_factor = factor_data.iloc[:, :]
y_factor = factor_data.index
# Using Kaiser-Meyer-Olkin (KMO) value to validate adequate sample size
kmo_all, kmo_model = calculate_kmo(X_factor)
print(format(kmo_model,".4f"))
In order to conduct exploratory factor analysis, the Kaiser-Meyer-Olkin (KMO) value must be greater than 0.5, as this shows that the sample size is fit for factor analysis. The closer the value is to 1, the better. Although the below calculated KMO value isn't ideal, it indicates that factor analysis is viable. This value could be improved upon with a bigger sample size, i.e. by collecting more data.
An eigenvalue returned by factor analysis that is greater than 1, is considered as selection criteria for the feature. This is because the eigenvalue represents the amount of variance created by that variable. So an eigenvalue with a value of 4.04, for example, accounts for variance created by 4.04 variables.
# Create scree plot using matplotlib
plt.scatter(range(1, X_factor.shape[1] + 1), ev)
plt.plot(range(1, X_factor.shape[1] + 1), ev)
plt.title("Scree Plot")
plt.xlabel("Factors")
plt.ylabel("Eigenvalue")
plt.grid()
plt.show()
The newly created latent factors, i.e. features created from the observable features:¶
Factor_1 = []
Factor_2 = []
Factor_3 = []
for idx, row in fa_loadings.iterrows():
if row.Fact1_loadings_ > 0.1:
Factor_1.append(row.Original_Factors)
elif row.Fact2_loadings_ > 0.1:
Factor_2.append(row.Original_Factors)
elif row.Fact3_loadings_ > 0.1:
Factor_3.append(row.Original_Factors)
print("Factor_1: ", Factor_1)
print("\nFactor_2: ", Factor_2)
print("\nFactor_3: ", Factor_3)
How can the viability of the returned factors be measured? With the Cronbach alpha, factors can be measured to determine whether or not the variables of a factor form a “coherent” and reliable factor. A value above 0.6 for the alpha is in practice deemed acceptable. Here the Cronbach alpha is calculated using the pingouin package:
# Get cronbach alpha
factor1_alpha = pg.cronbach_alpha(factor1)
factor2_alpha = pg.cronbach_alpha(factor2)
factor3_alpha = pg.cronbach_alpha(factor3)
print(
"factor1_alpha: ",
factor1_alpha,
"\nfactor2_alpha: ",
factor2_alpha,
"\nfactor3_alpha: ",
factor3_alpha,
)
Logistic Regression¶
Correlation Amongst All Skills¶
fig, ax = plt.subplots(figsize=(16, 16))
data_corr = data.corr()
cmap = sns.mpl_palette("Set3", 4)
# Generate a Heatmap
sns.heatmap(data_corr, cbar_kws={"shrink": 0.8}, square=True, cmap=cmap)
# title
title = "Correlation Between All Skills"
plt.title(title, loc="center", fontsize=24)
Visualizing Class Balance¶
It's important that classes be balanced because with strong class imbalance, a model can seem to perform well simply by predicting that every class is the majority class. But this is not indicative of overall model accuracy.
# Looking at Class Imbalance
sns.countplot(x="Frontend", data=data, palette="hls")
print(
"Percentage of frontend developers is: ",
"{:.2f}".format(len(data[data["profile"] == "Frontend"]) / len(data) * 100),
"%",
)
print(
"Percentage of non-frontend developers: ",
"{:.2f}".format(
(len(data) - len(data[data["profile"] == "Frontend"])) / len(data) * 100
),
"%",
)
plt.show()
Fitting the Logistic Regression Model to our Data¶
from sklearn.linear_model import LogisticRegression
logreg = LogisticRegression()
logreg.fit(X_train, y_train)
Measuring Model Performance¶
- Accuracy
- The Confusion Matrix
- The Classification Report
- ROC_AUC_score
- ROC Curve
Scoring Metrics¶
Accuracy can be a useful measure when the same amount of samples exist per class in a dataset. But with an imbalanced set of samples, accuracy is biased to the size of the test data. A more robust scoring metric is the f1 metric, used along with cross_val_score to score the model across multiple test runs.
# Generating predictions by which to measure accuracy of the fitted model
y_pred = logreg.predict(X_test)
print(
"The f1 score of the logistic regression classifier on test set: ",
"{:.3f}".format(
np.mean(cross_val_score(logreg, X_test, y_test, scoring="f1", cv=5))
),
)
The Confusion Matrix¶
This is the basis from which all other metrics are calculated. It provides a simple, but clear, picture which shows how the model is performing.
# Confusion Matrix
cm = metrics.confusion_matrix(y_test, y_pred)
group_names = ["True Negative", "False Positive", "False Negative", "True Positive"]
group_counts = ["{0:0.0f}".format(value) for value in cm.flatten()]
group_percentages = ["{0:.2%}".format(value) for value in cm.flatten() / np.sum(cm)]
labels = [
f"{v1}\n{v2}\n{v3}"
for v1, v2, v3 in zip(group_names, group_counts, group_percentages)
]
labels = np.asarray(labels).reshape(2, 2)
sns.heatmap(cm, annot=labels, fmt="")
The ROC Curve¶
This is a commonly used curve that summarizes the performance of a classifier over all possible thresholds. It's generated by plotting the True Positive Rate (y-axis) against the False Positive Rate (x-axis) as the threshold for assigning observations to a given class is varied.
plt.subplots(1, figsize=(10, 10))
plt.title("Receiver Operating Characteristic (ROC) : Logistic regression")
plt.plot(false_positive_rate, true_positive_rate)
plt.plot([0, 1], ls="--")
plt.plot([0, 0], [1, 0], c=".7"), plt.plot([1, 1], c=".7")
plt.ylabel("True Positive Rate")
plt.xlabel("False Positive Rate")
plt.show()
ROC_AUC Score¶
In order to return one number that reveals how good the ROC_AUC curve is, calculate the Area Under the ROC Curve, or the ROC_AUC score. This tells us how good the model is at ranking predictions.
# The Area Under the Curve (AUC) is the measure of the ability of a classifier to distinguish between classes
print(
"The roc_auc_score for the Logistic Regression classifier: ",
"{:.2f}".format(roc_auc_score(y_test, y_score)),
)
Logistic Regression with PySpark¶
sdf = sqlCtx.createDataFrame(data)
Transforming the Data to Spark¶
# Build a dict (orginal column name: new column name)
mapping = {col: col.replace(".", "_").replace(" ", "_") for col in selected_columns}
# Select all columns and create an alias if there is a mapping for this column
df_renamed = sdf.select(
[F.col("`" + c + "`").alias(mapping.get(c, c)) for c in sdf.columns]
)
# Create a VectorAssembler obj using the renamed columns as input
assembler = VectorAssembler(inputCols=list(mapping.values()), outputCol="features")
transformed_data = assembler.transform(df_renamed)
# Rename "Fullstack" column to "label"
data = transformed_data.withColumnRenamed("Fullstack", "label")
Imbalancing handling¶
In our dataset (train) we have 36.04 % negatives and 63.96 % positives. Since positives are in a majority. Therefore, logistic loss objective function should treat the negative class (Outcome == 0) with higher weight. For this purpose we calculate the BalancingRatio as follows:
Against every Outcome == 0, the calculated BalancingRatio is assigned to the corresponding row under the “classWeights” column. Likewise with every Outcome == 1, 1-BalancingRatio is assigned to the corresponding row in the “classWeights” column.
data = data.withColumn(
"classWeights",
F.when(data.label == 0, BalancingRatio).otherwise(1 - BalancingRatio),
)
Create PySpark Logistic Regression Model¶
from pyspark.ml.classification import LogisticRegression
# Create a logistic regression object
lr = LogisticRegression(
featuresCol="features", labelCol="label", weightCol="classWeights", maxIter=100
)
lrModel = lr.fit(train)
trainingSummary = lrModel.summary
Visualize Model Performance¶
# Plot the recall-precision curve
pr = trainingSummary.pr.toPandas()
plt.plot(pr["precision"], pr["recall"])
plt.xlabel("Precision / 1 - FPR")
plt.ylabel("Recall / TPR")
plt.show()
# Plot the threshold-precision curve
tp = trainingSummary.precisionByThreshold.toPandas()
plt.plot(tp["threshold"], tp["precision"])
plt.xlabel("Threshold")
plt.ylabel("Precision")
plt.show()
Binary Decision Tree¶
# Generating the DecisionTree model
classifier = DecisionTreeClassifier()
# Fitting the model
classifier.fit(X_train, y_train)
The Classification Report¶
While the below generated classification report reveals this Binary Decision Tree model to score an accuracy rate of 82%, a Multinomial (or Multi-class) model trained with exactly the same data and with the same model specifications will score an accuracy rate of 68% in predicting the five (5) classes present.
# Making predictions with the model
y_pred = classifier.predict(X_test)
# Generating the classification report
print(classification_report(y_test, y_pred))
The ROC Curve¶
from yellowbrick.classifier import ROCAUC
# Instantiating the visualizer with the Decision Tree model
visualizer = ROCAUC(classifier, classes=["Frontend", "Non-Frontend"])
visualizer.fit(X_train, y_train) # Fit the training data to the visualizer
visualizer.score(X_test, y_test) # Evaluate the model on the test data
visualizer.show()
Visualizing the Decision Tree's Strcuture¶
The tree will stabilize at a max_depth of 125 on its own, but for the sake of this visual, the max_depth has been set to 25.
clf = DecisionTreeClassifier(
max_depth=25
) # max_depth is maximum number of levels in the tree
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(14, 14), dpi=600)
clf = clf.fit(X, y)
plot_tree(clf, filled=True, feature_names=feature_names, class_names=class_names)
fig.set_size_inches(3, 6)
Multinomial Random Forest¶
Random Forests see high performance with data of high dimensionality. Furthermore, the algorithm handles outliers and unbalanced data well.
Fit the Random Forest Model to our Data¶
# Generate the RandomForest model
forest = RandomForestClassifier(class_weight="balanced", random_state=0)
# Fit the model
forest.fit(X_train, y_train)
The Classification Report¶
- Precision — What percent of your predictions were correct?
- Recall — What percent of the positive cases did you catch?
- F1 score — What percent of positive predictions were correct?
- Support is the number of actual occurrences of the class in the specified dataset.
from yellowbrick.classifier import ClassificationReport
visualizer = ClassificationReport(forest, classes=classes, support=True)
visualizer.fit(X_train, y_train) # Fit the visualizer and the model
visualizer.score(X_test, y_test) # Evaluate the model on the test data
visualizer.show()
The Multi-Class ROC Curve¶
from yellowbrick.classifier import ROCAUC
# Instantiating the visualizer with the Random Forest model
visualizer = ROCAUC(forest, classes=classes)
visualizer.fit(X_train, y_train) # Fit the training data to the visualizer
visualizer.score(X_test, y_test) # Evaluate the model on the test data
visualizer.show()
Conclusion¶
The original data's format required some creativite re-formatting in order to see it correctly processed by our Machine Learning algorithms. But ultimately, we were able to implement a binary classification problem with Logistic Regression & generate viable results with Multi-Class Random Forest model. Although these algorithms are powerful in their ability to make sense of data, it’s clear that the creativity and insight which only people can provide, will continue to be fundamental to the process of Machine Learning.
Future Work¶
In the future, I hope to investigate the information supplied by the predict.proba method in more depth. This method shows the probability that an observation belongs to a certain class. It sheds light on which samples the algorithm may be classifying with less certainty versus those classified with more certainty. A closer look at the skills of those profiles may lend insight as to which skills are determinative in predicting a profile type.